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Abstract 

In this paper we consider the Modified Craig-Sneyd (MCS) scheme which forms a promi¬ 
nent time stepping method of the Alternating Direction Implicit type for multidimensional 
time-dependent convection-diffusion equations with mixed spatial derivative terms. When the 
initial function is nonsmooth, which is often the case for example in financial mathematics, 
application of the MCS scheme can lead to spurious erratic behaviour of the numerical ap¬ 
proximations. We prove that this undesirable feature can be resolved by replacing the very 
first MCS timesteps by several (sub)steps of the implicit Euler scheme. This technique is 
often called Rannacher time stepping. We derive a useful convergence bound for the MCS 
scheme combined with Rannacher time stepping when it is applied to a model two-dimensional 
convection-diffusion equation with mixed-derivative term and with Dirac-delta initial data. 
Ample numerical experiments are provided that show the sharpness of our obtained error 
bound. 

Key words: Convection-diffusion equations, ADI splitting schemes, convergence analysis, Rannacher 
time stepping. 


1 Introduction 


In financial mathematics, the fair value u{si, S 2 ,t) of a European style option on two underlying 
assets is modelled by the two-dimensional Black-Scholes partial differential equation (PDE), see 

e g- n, 

Ut = \als\us^si + PCri(T2SlS2MsiS2 + \o-lslus^s2 + T’SlMsi + rS2Us2 “ (1-1) 


for Si, S 2 > 0, 0 < t < T. Here, t denotes the time to maturity T and we assume real parameters 
r, (Ti > 0, {72 > 0, IpI < 1. The PDE (1.1) is provided with an initial condition that is defined 
through the payoff of the option. 

The mixed spatial derivative term in (1.11 represents the correlation between both asset prices 
in the two-dimensional Black-Scholes model. Mixed spatial derivative terms are very important, 
notably, in the held of hnancial option valuation theory. Here they arise due to the correlation 
between underlying stochastic processes. 

A well-known approach for determining the fair values u(si, S 2 , T) consists of numerically solv¬ 
ing PDE (1.1) by the method-of-lines, whereby one hrst discretizes in space and subsequently in 
time. In this paper we consider a uniform Cartesian grid and second-order central hnite difference 
schemes in space. This semidiscretization is second-order convergent with respect to the spatial 
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mesh width if the initial and boundary data is smooth, see e.g. [10]. When the PDE is multidi¬ 
mensional, then the application of classical implicit time discretization methods to the obtained 
semidiscrete systems can be computationally very intensive. In view of this, for the effective time 
discretization, operator splitting schemes of the Alternating Direction Implicit (ADI) type are 
widely considered. In this paper we consider the Modified Craig-Sneyd (MCS) scheme |S], which 
is a prominent scheme of the ADI type. In the past years various positive stability results for 
the MCS scheme have been derived relevant to multidimensional convection-diffusion equations 
with mixed derivative terms, see e.g. HHIIHIIIII. Recently, in ’t Hout and Wyns [5] proved 
that, under some natural stability and smoothness assumptions, the MCS scheme is second-order 
convergent with respect to the time step whenever it is applied to semidiscrete two-dimensional 
convection-diffusion equations with mixed derivative term. The temporal convergence result from 
[9| has the crucial property that it holds uniformly in the spatial mesh width. Hence, the fully 
discrete numerical solution is second-order convergent in space and time for smooth initial and 
boundary data. 

A relevant convergence analysis for the MCS scheme and nonsmooth data is still open in the 
literature. In financial applications, however, the initial function is in general nonsmooth. It 
is well-known that convergence can then be seriously impaired. As an illustration, consider a 
two-asset cash-or-nothing option with strikes ATi > 0 and K 2 > 0, so that 


m(si,S 2,0) = l{si>Ari}l{s2>A^2}i 

where 1 denotes the indicator function. In the upper left plot in Figure the numerical solution 
for m(si,S2,T) is shown for (natural) financial parameter values r = 0.05, cri = 0.2, (T2 = 0.25, 
p = —0.7, Ki = \, K 2 = I, T = 2. Irregularities can be observed around the strikes, leading to a 
loss of accuracy in the maximum norm. For hedging purposes it is important to consider also the 
Greeks, for example the cross gamma V = Usiss- The corresponding PDF is given by 

Ti — 2^1 ^iTsiSi T s\S2 A 2^2^2^S2S2 

+ {r + af+ paia 2 )siTs^ + {r + + pcria 2 )s 2 Ts 2 + (r -I- paia 2 )T, (1.2) 

for si, S 2 > 0, 0 < t < T. This is supplemented with initial function 

r(si,S2,0) = UsiS2(si,S 2,0) = 5(si - Ki)S{s2 - K 2 ), 

where 6 is the Dirac delta function. The lower left plot in Figure[2shows the numerical solution for 
the cross gamma at maturity T for the same financial parameter values as above. Around the point 
(si,S 2 ) = (All, 7 ^ 2 ) strong, spurious erratic behaviour shows up and, hence, this approximation 
is useless in practice. If the cross gamma is approximated by applying finite difference schemes 
directly to the numerical solution for the option value, which is a common alternative technique 
in practice, the same observations are found. 

For one-dimensional applications in finance, the impact of nonsmooth initial data on conver¬ 
gence has already been studied extensively and various techniques have been proposed in order 
to recover standard convergence results, see e.g. laEU. A common technique consists of first 
applying several implicit Fuler (sub)steps and then continue with the time stepping scheme under 
consideration, [I3j . This is called Rannacher time stepping or implicit Euler damping. 

Consider again PDFs (0) and 0 for the two-asset cash-or-nothing option. Replacing the 
MCS scheme in the first two timesteps by four half-timesteps of the implicit Fuler scheme, the 
two right plots in Figure]^ are obtained. Clearly, there are no longer irregularities or oscillations 
present. In many other multidimensional applications, see e.g. @|, the same observations were 
made. To the best of our knowledge, however, there are no theoretical results available in the 
literature concerning the favourable effect of Rannacher time stepping on the convergence of the 
MCS scheme if the initial data is nonsmooth. 

In the present paper we will prove a useful convergence bound for the MCS scheme when it is 
applied to a model two-dimensional convection-diffusion equation with mixed derivative term, pro¬ 
vided with Dirac delta initial data. Here, semidiscretization is performed with second-order central 
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No Rannacher time stepping Rannacher time stepping with four haif-timesteps 



Figure 1: Numerical approximations of the cash-or-nothing option value (top) and of its cross 
gamma (bottom) without (left) and with (right) Rannacher time stepping with four half-timesteps. 
The financial parameter values are r = 0.05, ai = 0.2, cr 2 = 0.25, p = —0.7, Ki = 1, K 2 = 
1, T = 2. 


finite difference schemes. The precise influence of Rannacher time stepping on the order of conver¬ 
gence will be investigated. Our analysis in this paper is inspired by that of Giles and Carter [3], 
who deal with the Crank-Nicolson scheme applied to a model one-dimensional convection-diffusion 
equation. We make use of a two-dimensional mixed discrete/continuous Fourier transformation 
and analyse the asymptotic behaviour of the Fourier transform. By applying then the inverse 
transformation we arrive at an error bound for the total error. The sharpness of the error bound 
is confirmed by ample numerical experiments. 


3 




























2 The Modified Craig—Sneyd scheme 


Semidiscretization by finite difference methods of initial-boundary value problems for time-dependent 
convection-diffusion equations leads to large systems of stiff ordinary differential equations (ODEs), 

U'{t)=F{t,U{t)) (0<t<T), U{Q) = Uo, 

with given operator F and given initial value Uq. Assume the PDE is two-dimensional and the 
semidiscrete operator F is decomposed into a sum 

F{t,v) = Fo{t,v) + Fi{t,v) + F 2 {t,v) (0 < t < T), 


where Fq represents the mixed spatial derivative term and Fi, F 2 , represent all spatial derivative 
terms in the first, respectively, the second spatial direction. Let 0 > 0 be a given parameter, 
> 1 the number of timesteps and set = nAt with At = T/N. Then the Modified Craig- 
Sneyd (MCS) scheme generates, in a one-step fashion, approximations Un to U{tn) successively 
for n = 1, 2 ,..., A through 


Yo — Un-l + At F{tn-l, Un-l), 

Yi = Yi^i + 9At{F,{tn,Yi) - Fi{tn-i, C/„-i)), * = 1,2, 

lo = Iq + (Foitn, Y 2 ) — Fo(t„_i, Un-l)) , 

< .. ^ (2-1) 
^0 = lb + (5 - 0)At {F{tn, Y 2 ) - F{tn-l, Un-l)) , 

Yi = Yi-I + eAt{Fi{tn,Yi) - Fi{tn-1, Un-l)), * = 1,2, 

. Un=%- 


The MCS scheme (2.1) was introduced by in ’t Hout & Welfert [8] for general multidimensional 
convection-diffusion problems with mixed derivative terms. It can be viewed as an extension of 
the Craig-Sneyd (CS) scheme, proposed in [5]. For 9 = 1/2, the MCS scheme reduces to the CS 
scheme. Besides 9 = 1/2, common choices for 9 in the literature are 0 = 1/3 and 0 = 1. Scheme 
(2.1) starts with an explicit Euler stage applied to the full system, which is followed by two implicit 
corrections corresponding to each of the two spatial directions. Subsequently an explicit update 
is performed, followed again by two implicit unidirectional corrector stages. Note that both F 
and Fq, which contain the mixed derivative term, are always treated in an explicit manner. Each 
implicit stage handles spatial derivatives in only one spatial direction. This can lead to a major 
computational advantage in comparison to classical non-splitted implicit time stepping methods. 


3 Model problem 


Consider the coordinate transformation x = •\/21og(si)/(Ti and y = •\/21og(s2)/cr2- The PDE (1.11 
is then transformed into 


Ut — Uxx + 2pUxy Uyy + ('^ 

for —00 < Xi, a ;2 < 00 , 0 < t < T. This provides a motivation for considering a constant coefficient 
model convection-diffusion equation with mixed derivative term 


V2r 


Ut - Uxx F *^pUxy F Uyy -(- CLlUx -(- U2Uy, 


(3.1) 


for —00 < x,y < 00 , 0 < t < T = 1 and with \p\ < 1. We supplement equation (3.1) with the 
initial condition 

u{x,y,Qi) = S{x)S{y), 
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which arises for example in the case of the cross gamma of a two-asset cash-or-nothing option. 
The Dirac delta initial function, however, has other important applications as well. For instance, 
it arises naturally in the adjoint equation for the joint density. By using the Fourier transform 
pair 


/ OO nOO 

/ u{x, y, t) exp(— Iko;) exp{— ir]y)dxdy, 

-OO j —OO 
^ pOO nOQ 

u{x,y,t) = / / u{K,r],t) exp{iKx) exp{iriy)dKdri, 

J — OO J — OO 

an exact closed-form analytical solution will be derived. Here i denotes the imaginary unit. Taking 
the Fourier transformation of equation (3.1) yields the ODE 

Ut = —K^u — 2pKT]u — rfii + ia\KU + ia2riu, 

subject to initial condition u{k, 77 , 0) = 1. The solution of this transformed equation is given by 

u{k, rj, t) = exp(—-|- 2pKri + rj^ — ioiK — ia 2 T])t). (3-2) 

Next, if {Xi,X 2 ) is a multivariate normal distributed random variable with mean (/ri,/i 2 ) and 
covariance matrix E, its characteristic function is defined by 

E[exp(iKXi) exp{ir]X2)] = exp(iK/ri -|- ir]pL2 — y)^)- 

By exploring the connection between the characteristic function of a random variable and the 
Fourier transform of its density function, it follows that u{x, y, t) can be seen as the density function 
of a two-dimensional normal distributed random variable with mean (/ii, 7 x 2 ) and covariance matrix 
S given by 

{di,d 2 ) = {-ait,-a 2 t) and Y. = 

Since \p\ < 1, this yields the closed-form analytical solution 

u{x,y,t) = ^-^-^^exp(^-j^j^[{x + ait)‘^ + {y + a 2 t)‘^ -2p{x + ait){y + a 2 t)]^ . 

4 Discretization 

As mentioned in Section spatial discretisation of (3.1) will be performed on a uniform Carte¬ 
sian grid with second-order central finite difference schemes. For the time integration the MCS 
scheme will be considered. Let hi denote the spatial mesh width in the a;-direction, /12 the spatial 
mesh width in the y-direction and define spatial gridpoints {xj,yk) = {jhi,kh 2 ) for all j,k € Z. 
Semidiscretization of ( |3.1| ) with second-order central finite difference schemes then gives rise to 
approximations Uj^k(t) of the exact solution value u{xj,yk,t) which are defined by the system 


= AC/,. fc(t). 


(4.1) 


where A = Aq -I- Ai -|- A 2 and 


Aq 

Ai 

A 2 


2hih 


-^2x^2y, 


1 c2 , A 


02 
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with S 2 X, Sy the usual second-order central finite difference operators. For example, 

^xUj,k{t) = Uj-i^k{t) — ‘2Uj^k{t) + Uj + i^k{t), 

^2x^2yUj^k(i) — + (^) “t“ — — l —(0 ■ 


Semidiscrete system (4.1) is provided with initial data 

UjAo) 


hih2 

0 


if j = fc = 0 , 

else, 


in order to approximate the Dirac delta initial function. For convenience we define 

Z = AtA, Zi = AtAi for i = 0,1,2, 

and we denote by I the identity operator. Then, starting from t/oj,fe = Uj^k{0), application of the 
MCS scheme to semidiscrete system (4.1) yields approximations Un,j,k of Uj^k{tn) successively for 
n = 1,2,... ,N through 

^0,j,k — d” ^')Un—l,j,ki 

{I — 9Zi)Yij^k = Yi-lJ.k — dZiUn-l,j,k 

boj',fc — Z^ j k 9ZQY2,j,k 9ZQUn—l,j,k-; 

%,,,k = %,j,k + {.\-9)ZY2,j,k-{\-e)ZUn-l,j,k, 


f = 1,2, 


(4.2) 


* = 1 , 2 , 


Un,j,k — F2,j,fc- 

Concerning the Rannacher time stepping, let Nq denote the number of initial MCS time steps 
replaced by 2Nq half-time steps of implicit Euler integration. Whenever Nq > 0 scheme (4.2) is 
replaced by 

2^')^n—ll2,j,k = Un—l,j,ki 

(4.3) 

for n = 1,2,..., min{7Vo, A*}. This provides a numerical approximation C/jv of the exact solution. 
The goal of our convergence analysis consists of quantifying the total error 

UN,j,k - u{xj,yk,l). (4.4) 

To do so, we will analyse the asymptotic behaviour of a mixed discrete/continuous Fourier trans¬ 
form for hi,h 2 ,At simultaneously tending to zero. Applying the inverse Fourier transformation 
on the resulting error in Fourier space will yield a useful bound for the total error (4.4). Special 
attention will be paid to the influence of Nq, i.e. the influence of Rannacher time stepping, on the 
total error. 


5 Asymptotic analysis in Fourier space 


We consider a mixed discrete/continuous Fourier transform pair, cf. e.g. 


00 00 


F(i9i,i?2) = A 1/12 X! ^i.fcexp(-ijrli)exp(-ifci?2), -tt < di,■d2 < tt, 

j — — oo k— — oo 


Vj,k = 


1 


/*7r pTT 


47r2/li/l2 


V{'&i,'92) exp(iji5i) exp(ifci?2)d*?id*^2. 


j, k gZ. 
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For ease of presentation, the dependency of the Fourier transform on 'di and '^2 will be omitted 
in the notation. 

Fourier transformation of t/oj'.fc yields Uq = 1. Concerning operator Zq it follows that 


OO OO 


ZqV = hih2 ^ ZoVj^kexp{-ij'di)exp{-ik'd2) 

j — — oo k— — oo 


OO OO 

= ^ Y (fj+i.fe +1 + - Fj-i,fc+i) exp(-ijdi)exp(-i/cd 2 ) 

j — — oo k— — oo 


OO OO 


= ^exp(idi)exp(ii? 2 ) Y X! ^i+i.fc+i exp(-i(j + l)di) exp(-i(fc + 1 )^ 2 ) 

j = — OQ k — — <X) 


OO OO 


+ ^exp(-idi)exp(-id 2 ) Y X! exp(-i(j - l)i9i) exp(-i(fc - l)'d 2 ) 

j——OQ k— — (yci 


00 00 


- ^exp(idi)exp(-id 2 ) Y X! Kj+Ffc-i exp(-i(j + l)di) exp(-i(A: - 1 )^ 2 ) 

j— — OQ k— — (yci 


00 00 


- ^exp(-ii?i)exp(id 2 ) Y X! Kj-i.fe+i exp(-i(j - l)di) exp(-i(A: + 1 )^ 2 ) 

j— — OQ k— — (yci 

= [exp(idi) exp(i'd 2 ) + exp(-idi) exp(-id 2 ) - exp(idi) exp(-i'd 2 ) - exp(-idi) exp(id 2 )] V 

= -l?f(sindisind2)F. 

Analogously one finds 


ZiF = sin^ ^ + iaif^ sindi^ V, 


sin2 ^ + ia2#^ sind2 ) V. 


Define functions 


(- 

zo = zo(i9i,i92) = sindi sin'd 2 , 

zi = zi(di,'d2) = -^sin^ + ioi^sindi, 

Z2 = ^2(di,'d2) = -^sin^ ^ + ia2^sind2, 


and z = zq + Zi + Z 2 . Then, Fourier transformation of the implicit Euler scheme (4.3) gives 

2 


I 7 r,= 


1 

1 - ^z 


Un-l- 


After some calculations, Fourier transformation of the MCS scheme (4.2) yields 

Un — 

with 


o 1 , ^ , (^^0 + (2 - 

H = I -\ -^--- -^5 - 

p 


where 


p = (1 - 6 > 2 :i )(1 - dz2). 

Assume that Nq < N. Since C/q = 1 it follows that 

2No 




(5.1) 


(5.2) 
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By applying the inverse Fourier transformation, the numerical approximation at t = T = 1 can 
be written as 


UN,j,k = 


1 


47r2/lift,2 J-nJ—. 


UN{'&i,'d 2 ) exp(ij-(9i) exp{ik'd 2 )d'did'd 2 


1 


^7r//i2 pT^lhi 


' —■Klh 2 —tt/Zii 

where we made use of the substitutions 


? 7 Ar(K/ii, 77/12) exp(ia;jK) exp{iykr])dKdr], 


di = nhi, 7/2 = 11^2- 

From Section it can be seen that the exact solution is given by 

1 


^00 pOO 


/ oo 

u{k, 77 , 1 ) exp(ia;«;) exp{iyr])dKdr]. 

-00 


In our analysis, we will examine the Fourier error 


/7Ar(/c/ii, 77 / 72 ) — m(k, 77 , 1 ) for — tt < 7 c/ 7 i, 77/12 < tt. 


(5.3) 


For hi,h 2 tending to zero, the total error (4.4) is approximated by 


2 r■^/h2 l-Tr/hi 

—^ / / (/7Ar(K/7i,77/72)-«(«:,77, 1)1 exp(ixjK)exp(i7/fe77)(i7t(i77. ( 5 . 4 ) 

47r J-Tr/h2 J-Tv/hi ^ ' 


Note that expression (5.4) can be viewed as the inverse mixed discrete/continuous Fourier trans¬ 


form of the Fourier error (5.3). 

In Figure juj is shown in the (di, 7?2)-domain for parameter values p = —0.7, Oi = 2, 02 = 3. 
This has to be compared with Figure where \Un\ is shown for the same parameter values. 
Discretization is performed with hi = h 2 = 1/6, At = 1/8 and well-known MCS parameters 
9 = 1/3,1/2,1. For the Rannacher time stepping we considered values Nq = 0, 2. From Figure]^ 
and Figure it is clear that the difference Un — u has different properties in different regions of 
the Fourier domain. These regions are illustrated in Figure]^ 

First there is a low-wavenumber region ®, where both |7?i| and |7?2| are small, in which there 
is a good agreement between Un and u. Next, if either |7/i| or I 7 / 2 I is medium and the other one 
is small or medium (region (2)), then both the Fourier transforms of the numerical solution and 
analytical solution are negligible. In the high-wavenumber region (3), i.e. where both |7?i|, 17 / 2 ! are 
large, we observe that the modulus of the Fourier transform u is close to zero. The modulus |/7jv|, 
however, is strongly dependent on Nq and the MCS parameter 9. For larger values of 9 we see 
that Um has a larger magnitude in the high-wavenumber region. Hence, a larger high-wavenumber 
error can be expected for larger values of 9. Further we observe that the modulus of Un in the 
high-wavenumber region is always damped whenever Rannacher time stepping is applied. This 
matches our observations from Figure where unwanted erratic behaviour was avoided by using 
Rannacher time stepping. Finally, we have the case where either |7/i| or |7?2| is large but the other 
one is not. In our analysis, the region @ where |7?i| is large and the region (5) where |792| is large 
will be treated separately. In both regions the Fourier transform u is negligible but Un has to be 
further analysed. In particular, we will show that Un is not negligible if the MCS scheme reduces 
to the CS scheme. ^ 

Following Giles & Carter [3] we will perform an asymptotic analysis of the Fourier error Un — u 
in each of these (five) disjoint regions which form a partition of the Fourier domain. We consider 
the limit hi, h 2 ,At —>■ 0 and since the same discretization is performed in both spatial directions. 


c = / 72//71 
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Figure 2: Magnitude of the Fourier transform of the exact solution u{x^ y, 1) for parameter values 
p = —0.7, Oi = 2,02 = 3. 


is held fixed. For ease of presentation we denote h = hi. Further, since both the semidiscretization 
and the time integration are convergent of order two for smooth initial data, it seems natural to 
keep 

A = At/h 

constant. Substitutions di = Khi,‘d 2 = r]h 2 yield 


zq = —sin k/i sin cTy/i = —^(cos((k — cri)h) — cos((k + crfjh)), (5.5a) 

zi = — ^ sin^ ^ + ioiAsinK/i = — ^(1 — cos Kh) + ioiAsinK/i, (5.5b) 

Z2 = sin^ ^ + ia2^ smcrjh = —^(1 — coscrjh) + ia2^ sincr/Zi. (5.5c) 


The expressions in (5.5) will be used to analyse the asymptotic behaviour of (5.2) as —)■ 0. 
Throughout the analysis, by the notation O {f{K,ri,h)) we shall always mean that the modulus 
\ ■ \ of the term under consideration is bounded by a positive constant times f{K, rj, h) where the 
constant is independent of k, ij and the mesh width h. In order to deal with the powers in expression 
(5.2) a log-transformation of t/jv will be considered. Since T = 1, thus N = l/(A/i), it follows that 


logC/AT = (iV-Afo)log(i?) + 2iVolog(^Y3^) (5.6a) 

= ^ [log(p2 +pz + 9zoz + (i - 9)z^) - 2 log(p)] (5.6b) 

+ No [2 log(p) - log(p2 +pz + 9zoz + (i - 9)z^) - 2 log(l - ^z)] . (5.6c) 


5.1 Taylor expansion oi Un 

Multiple regions will encounter values |«;|, \cr]\ < h~'^ with certain q < 1/2. By Taylor expansion 
of (5.5a) it directly follows that 




pX f {k + cr])^h^ (k — cr])^h^ {k -|- cr])‘^h'^ (k — cr])^h‘^ 


4! 


4! 


+ ■ 


— — ^Kcrjh — ^(«^^ + (?rj^)ncphf’ ^ 


-I- z^o^^ T- z^o^h^, 
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Time marching: 0 = 0.333 , N 


Time marching: 0 = 0.333 , N 


Time marching: 0 = 0.5 , N 


Time marching: 0 = 0.5 , N 


Time marching: 0=1 , N = 0 


Time marching: 0 = 1 , N = 2 


Figure 3: Magnitude of the Fourier transform Un with A^o = 0 (left) and Nq = 2 (right) for MCS 
parameter 6=1/3 (top), 0 = 1/2 (middle) and 6=1 (bottom). The other parameter values are: 
p = —0.7, oi = 2,02 =3,hi = 0-2 = 1/6, At = 1/8. 








Figure 4: Illustration of the different disjoint regions of the Fourier domain. 


where 


'0 

d3] 

'0 


= —2pXKr], 

= + c^rf)Kri, 


3^ 

l^od < 


Analogously as above, Taylor expansion of (5.5b I and (5.5c) yields 

Zi{h) = h + , 

Z2{h) = zl^^ h + z^^'^ + Z2^ , 


where 


= —Xk^ + IoiAk, 

= ^Xk‘^ — ^'iqiXk^, 

< IfA/i® + A|adA|«:|3 


jsli 


J5]| 


—Ayy^ + \a2Xp, 
^Xc^p'^ — ^ia2Xc^p^, 


< ^Xc*p^ + ^\a2\Xc^ 


Since g < 1/2, it is ensured that all terms in the above expansions stay bounded as h tends to 
zero. Using these expansions and the definition (5.1) of p it follows that 


p(h) = l+pW/i + p[2]/i2_^p[3]/j3_^p[4]^4_^p[5];^5^ 
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where 


pW = _0(^W+zW), 

p[2] = 

pm = _0(43i+43]^^ 

pw = 0 ^(zWzr+.f'zW), 

pm = 0 (l + (^2 + 0 ^ 77 ^)^) . 


Under the condition |k|, \cri\ < h~'^ with certain q < 1/2, the variables k and q can become very 
large as h tends to zero. In this case the highest powers of k, 77 will dominate the order term in 
pm. Under the same condition, however, k and 77 can both be very small and then the lowest 
powers of n, 77 will dominate. By considering the sum of 1 and the highest powers of k, 77 in the 
remaining order term, we ensure that both cases are covered. 

As mentioned above we will make use of log-transformation (5.6a) to analyse the asymptotic 
behaviour. Let / be a strictly positive and sufficiently smooth function and set 


g{h) = log(/(h)) for h > 0 . 


Taylor expansion yields 


g{h) = log(/(0)) + gmh + gmh^ + + g^h^ (5.7) 


where 


7 W 


7[2] 


7[31 


.[4] 


/'(O) 
/(O) ’ 


K 


/"(o) _ fM 


/(O) /(O)- 


/ 


( f"'{0) _ q /'( 0 )/"( 0 ) , o/'O)"/ 

V /(o) f{or +^/(0)3/ 


y f _ 4/'(U/"'(U+3/"(U" 
4! V /(U /(O" 


1 o /'(UU"(U _ afiOl'] 

/(?)" ° /(U^ ) 


with certain 0 < ^ < h. 


In order to encounter the first part of (5.6b|) consider 


fM{h) = p{hf +p{h)z{h) + 9zo{h)z{h) -b (5 - d)z{hf, 


so that 

/mW = 2p{h)p'{h) -b p'{h)z{h) -b p{h)z\h) + 6»(zq(/7)z(/i) -b zo{h)z\h)) -b 2(i - e)z{h)z'{h), 

/m(^) = 2p'(^)^ + ‘2pih)p''{h) + p"{h)z{h) -b 2p’{h)z’{h) + p{h)z”{h) 

+ 9{zQ{h)z{h) + 2zQ{h)z'{h) + Zoih)z”{h)) + 2(i — 0)(z'(h)^ + z{h)z”{h)), 

ImW = ^P {h)p”{h) + 2p{h)p"{h) + p”{h)z{h) + 3p'{h)z'{h) + ip {h)z”{h) + p{h)z"’{h) 

-b 9{z'q {h)z{h) -b 3zQ{h)z'{h) + 3zg(h)z''(h) + zo{h)z'”{h)) 

+ 2(i — 9){3z'{h)z"{h) -b z{h)z'"{h)), 

/mH^) = O + , 

and thus 


/m(0) = 1, 

//^(O) = 2 p'( 0 ) + U( 0 ), 

ritiQ) = 2p'(O)2 + 2p"(O)+2p'(O)U(O) + 20z'(O)U(O) + 2(i-0)U(O)^ 
/)(^(0) = 6p'(0)p"(0) + 2p'"(0) + 3p"(0)z'(0) + z"'(0). 
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Concerning the Rannacher time stepping, define 


fNoih) = 1 - lz{h), 

such that 

fNoiO) = 1 , 

/W(0) = -i 0 «(O) for * = 1,2,3, 

C^{h) = 0{l + {K^ + cWy). 

All of these expressions will be used in the forthcoming subsections, where we analyse the asymp¬ 
totic behaviour of t/jv in five different regions of the Fourier domain, i.e. the (k, ? 7 )-domain with 
\k\, \t]\ < T^/h. 


5.2 Region 1: |k|, |c? 7 | < h with q < 1/3 

In order to analyse log Un in this region, the parts stemming from the MCS scheme and Rannacher 
time stepping will be considered separately. Write ( |5.6b[ ) as 

XS [\og{p^ +pz + Ozqz + (i - 6»)z2) - 21og(p)] = ^ [\og{fM{h)) - 21og(p(h))] . 

Using the analysis above it follows that 

^ [logifMih)) - 21og(p(/r))] = sM + slU/r + (5.8) 

where 

As[°l = 2p'(0)-hz'(O)-2p'(0), 

Asl^l = i [2p'(0)2 -h 2p"(0) + 2p'(0)z'(0) 26lz(,(0)z'(0) + 2(i - e)z'{0f - (2p'(0) -f z'(0))2] 

- [p"( 0 )-p'( 0 )^] 

As[ 2] = 1 [6p'(0)p"(0) -f 2p'"(0) -h 3p"(0)z'(0) + z"'{0) 

- 3(2p'(0) -f z'(0))(2p'(0)2 2p"(0) -h 2p'(0)z'(0) 26»4(0)z'(0) -h 2(i - e)z'{0f) 

+ 2(2p'(0) + z'(0))3] - i [p"'(0) - 3p'(0)p"(0) + 2p'(0)3] 

= o{l + {K^ + cWf)- 


By using the expansions in Subsection 5.1 and after simplifying the resulting expressions, one gets 


^[o] 

/i] 

,[ 2 ] 


= — 2ptvr] — iaiK + ia 2 ^ 7 , 

= 0 , 

= + Ip{k^ + c^p^)Kri + | | ia2C^p 

— X^9^{—k^ + iaiK){—ri^ + ia2p)(—K^ — 2pKp — rf + ioiK + ia2p) 

>? ( u-'Z ..,2 I • ^ I ■ 


+ - 2pKp - *7^ + iaiK -I- ia2?7)'^ 

— A^(—— 2pKp — + ioiK + ia2p)i—pKp + {\ — 0 ){—k^ — + \aiK + ia 2 ? 7 ))^ 


As for the part stemming from the Rannacher time stepping, write (5.6c) as 


A^o[21og(p) - \og{fM{h)) - 21og(/jvo(/i))]. 


Using the same analysis as above one gets 

iVo[21og(p) - \ogUM{h)) - 21og(/jvo(h))] = (5.9) 
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where 


N: 


[ 1 ] 


4^’ 


No 

-As[°l -2(-iz'(0)) 

= 0, 

No 

-AsW-(-iz"(0)- 

iz'ior) 


N, 


[3] 


= \NoX‘^{-k‘^- 2pKr]-ri‘^ + iaiK+ia2r])‘^, 

= O {I + {k'^ + c^r,'^)^) . 


By combining (5.6a), (5.8) and (|5.9[) it directly follows that 


\ogUN = — 2pKr] — + iaiK + ia2P + + iVg )/i^ + + iYg )/i^ 


and hence 

Un = exp(—— 2pKp — ? 7 ^ + ioiK + ia 2 p) exp((sf^] + N'^^)h'^ + (s^^l + N'^^)h?). 


Next, we will expand the second exponential in order to compare this expression with the Fourier 
transform u from (3.2) at t = 1. Let 

e{h) = exp(c^^^h^ + 


where 

then 


cl^l = 0 (^1 + (^2 c 2 ^ 2 ^ 3 ^ ^ ^[ 3 ] ^ (^2 , 


e'(h) = (2c[^]/i + 3cf^]/i^)e(/i), 

e''{h) = (2c[^] + 6cf^]/i)e(/i) + (2cf^]/i + 3cf^]/i^)^e(/i), 

e'''{h) = &Se{h) + 3(2cl^' + Qc'^^h){2c^'^^h + ■iSh'^)e{h) + {2c^'^^h + -iSh^fe^h). 


Since |k|, |cry| < h with q < 1/3, we have that e(0) = 1 and \e{h)\ < exp(l) whenever h is 
sufficiently small. Hence it follows that 

e{h) = 1 + 

with 

el^l = cl^l, 

e[3l = C>(1 + (k2 + c27?2)4) + 

= O (1 + + cW)^) + O (1 + (ac^ + cWf) h, 

where the latter equality follows from the assumption |k|, \cp\ < h~‘^ with g < 1/3. Finally, for 
this region, one arrives at the following expression for the Fourier error ( |5.3[ ): 

h^u{K, p, 1) ((st^l + Tvf') + 0{1 + {k^+ C^p^f) h + 0 (l + {K^ + c^p^f) . (5.10) 

Note that and actually depend on k and p. For ease of presentation, this is omitted in 
the notation. 


5.3 Region 2: |k| < h |cr/| < h with qi,q 2 < 1/2 and with qi > 1/3 
or q 2 > 1/3 


First consider the case where both qi < 1/2 and q 2 < 1/2. 
5.2 expression (|5.6b ) can be rewritten as 


Based on the analysis in Subsection 


Nlog{R) = ^ [\og{p^+pz + ezoz + (i - e)z^) - 21og(p)] = s[°l + 
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where 


5 ^°^ = —K^ — 2 pKri — rf + laiK + 10277 , 

s[ 2 '] = 0{{K^+cWf)- 

Since either k or 77 becomes large in this region as h tends to zero, only the highest powers of k, 77 
are taken into account in the order term in 1. From \p\ < 1 one gets 

+ 2 pKr] + 77 ^ = (1 - | p |)( k ^ + 77 ^) + \p\{k + sgTi{p)r]f > 0 , 

such that 

7?.(s[°l) < -(1 - \p\){k^ + rf) < 0. 

Using that both both <71 < 1/2 and < 1/2, it directly follows that 

lim(K^ + c^rf‘)‘^h‘^ = 0 , 

7i->-0 

and thus 

3(5>0 3/io>0 W h < Hq : TZ{N\og{R)) < —S{k^ + if). 

Hence, for h < ho 

and since |k| > or |c 77 | > h~^l^ we may conclude 

\R^\ = 0{h^) Vw>0. (5.11) 

Next, consider the case where at least one of the equalities, qi = 1/2 or 52 = 1/2, holds. For 
analysing the asymptotic behaviour of R we then make use of the following proposition. Its proof 
is a direct modification of the proof of one of the statements in [3 Theorem 1] and is therefore 
omitted. 

Proposition 5.1 Let zo,zi,Z 2 denote real numbers with 

Zl<0, Z2<0, |3()| < 2\p\y/ziZ2, (5.12) 

and \p\ < 1. Set z := Zg + zi + £2 and p := {1 — 0zi)(l — Hi/)- /f U < 0 or Z 2 < 0, then 

^ +p£+ezo£ + {\ - 


whenever 0 >\ and 0 > 

Recall that in the current region of the Fourier domain the assumption |«;| < , \cr]\ < 

with < 71,52 < 1/2 holds. This yields 


lim Zglh) = lim —2p\Krih =: zn G K, 
Ti-J-O /i-J-O 

lim zi(h) = lim —Xn^h =: z) € 

h^O Ti-s-O 

lim Z 2 ih) = lim —Xp^h =: Z 2 € R~- 

Ti-fO Ti-s-O 


Since |k| = 
Zg,Zi,Z 2 in 


h or |c? 7 | = h it follows that z* < 0 or zj < 0. Hence, all the assumptions on 
Proposition |5 .1 1 are fulfilled such that 


lim \R\ = 

Ti-J-O 






p^ +pz + OzqZ + (^ - d)z 


jp 


< 1 , 
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and thus 

M 

for 

Further, it always holds that TZ{z) < 0 such that 



(5.13) 

9> 1 and 9> . 

(5.14) 




By combining this with (5.11) and (5.13), and by using that A^o is independent from h, one may 
conclude that in this region it holds that 

|[7a,| = \R^\\R-^°\\1- =0{h'^) y-w > 0, 


under restriction (5.14) on 9. This means that |17Ar| quickly becomes negligible as h tends to zero. 
It decays faster to zero than any polynomial in h. 


5.4 Region 3: |k|, |c? 7 | > h with q > 1/2 

Here we reconsider the substitutions = nhi = nh, '^2 = ^ 7^2 = rjch in order to get 


Zo = 

-2p^ sindi sin ^ 2 , 

Zi = 

—4^ sin^ ^ + ioiAsindi, 

Z2 = 

-4^ sin^ ^ + 102^ sind2 


Further, in this region i9i, ^2 are different from zero and since we consider values —tt < ^2 < tt, 

we may write 


and thus 


16^ 


7? sin^ ^ sin^ ^ 


16 ^ sin^ ^ sin^ ^ 


16^ sin^ ^ sin" ^ 


2 ^ 


4^ sin^ % 


c cot ^ cot % 


zo = 


zi = 


2A 


-^2 = 


4A sin^ ^ 


4Asin2^ 


-(1 — 0 Z 2 ) = 9 + 


4Asin2^ 


4Asin2 ^ 


-h. 


cot % 2 

101 - ^TZi-h , 

SAsin^^ 

C cot ^ 2 

102 - r) q O , 

SAsin^^ 
cot ^ \ 


- 0ia2 


ccot ^ 


le^sin^ ^sin^ ^ 


p = 9^+9 


+ 


4Asin2 ^ 


cot ^ 


4Asin2^ 


- 6»ia2 


ccot ^ 


4A sin^ ’|l 


— dioi 


cot ^ 


4A sin^ 


- 0102 


ccot ^ 


h^. 
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Making use of an expansion similar to (5.7) it follows that 



le^sin^ ^sin^ ^ 


{p^ + pz + OzqZ + {^ -e)z^) 


4Asin2 ^ 


cot ^ 

- + 


„ c cot ^ cot ^ „ 

- e^p -1,— ^--e^ 


2 

1 


4Asin2 ^ 


2A 


4A sin^ ^ 


-e^ 


4Asin2 ^ 


ccot ^ 

-0102-^ 


h 

¥ 


+ O 



and 

log 


16^5 sin^ ^ sin^ ^ 


= log 02 


+ 


4Asin2 ^ 


cot ^ 

- eiai—¥- 


4A sin^ ^ 


- 0102 


ccot 



+ O 



Combining both expressions yields 

4A02 ' 


log {R) = - 


2pccot ^ cot ^ + 


+ O 



1 c2 sin^ ^ + 2pccos ^ sin ^ cos ^ sin ^ + sin^ ^ ^ 
sin2 ^ sin2 f ^ 



Further, recall that in this region i?i,d 2 are both different from zero such that 

1 ^ 2 ) := sin^ ^ + 2pccos ^ sin ^ cos ^ sin ^ + sin^ ^ > 0, 

and hence 


(5.15) 


logi?^ = NlogR 

_ 1 t(di,-02) 

4A20^ sin^ ^ sin^ ^ 

As for the implicit Euler time stepping scheme we note 

^ (1 — ^z) = (? sin^ ^ + ipcsini?! sind 2 + sin^ ^ — toi^ sindi — 102 ! sind 2 ^ h. 
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Using once again an expansion analogous to (5.7) it follows that 


sin^l 


log(l-iz) =log(|4)+log(6(l9i,^2))+ O 

which yields 

log = -27Vo log {§-) - 2No log ^ 2 )) + O 


I sin f I 


I sin ^ I 


I sin f I 


Making use of relationship (5.2) one becomes an expression for the logarithm of the Fourier 
transform Un ■ 


log Un = - 




sin^ ^ sin^ ^ 


1 + 0 



- 2iVolog (^) - 2Afolog(t(-di,-d2)) + O 


sm 


sm 


1?2 I 


such that in this region 


Un = 


[2Ai(i?i,i92)] 


2No 


exp 


1 


t(i?i,'d2) 


4A202 Sin2 sin^ % 


l + O 


(i^ii + i^2ir 


(5.16) 


In Figure]^ we noticed that in the high-wavenumber region, i.e. where both I'd!!, 11 ^ 2 1 a^re large, 
the norm \Un\ is highly dependent on the MCS parameter 9. This is confirmed by (5.16) since 
inequality ( |5.15 ) holds in the high-wavenumber region. Hence, for larger values of the MCS 
parameter 9 one can expect a larger high-wavenumber error. 

5.5 Region 4: |k| > |c? 7 | < with qi > 1/2, q 2 < 1/2 

Reconsider the substitution i?! = nh and recall that —tt < < tt. Then, as is non-zero in this 

region, one may write 


such that 


1 


l^sin^^ 


l^sin^^ 


■21 = 


■22 = 


1 


— 1 -|- iioihcot 

2 cr}h 


- „ . 2 2 

(? sm^ 


—R -h IO 2 


-—^^^sincTyh, 
4c sm 


1 

4^sin2^^ 


20 = 


— - cot ^ sincu/i, 
c 


= 0 


cot j /i j (1 - 9z2) , 


lim — 
h^o 4 ^ sin"^ ^ 

lim —;- 

h^o 4f sin^ ^ 

h 2 

1 


2 i9i 


22 = 0 , 


h^0 4^sin2 4^° 


lim 

x-s-O 

lim 
h^o 4 ^ sin^ 


h 2 

1 


0(1 -|- Z 2 ), 
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where denotes a positive real number. Hence, concerning R it follows that 


p^ +PZ + 9 zqz + {1 -9)z‘^ 


(1 - 1 - 22 )^ 0 ^ - (2 + Z2)9 + 4 

p2 


(1 + ^ 2)202 


lim 

h-^0 


For the latter expression we obtain the following positive result. 
Proposition 5.2 If 0 > 1/4 and 9 ^ 1/2, then 

1 


(1 + Z2Y0‘^ — (2 + Z2)9 + i 


< 1 

(1 + ^ 2)202 

for all real numbers Z 2 >0. If 9 = 1/4 or 9 — 1/2, then the inequality holds for numbers £2 > 0. 
Proof Let Z 2 be a positive real number. First, it is clear that the inequality holds whenever both 


— (2 + Z2)9 + ^ < 0, 

2(1 + 22)^^^ — (2 + ^2)9 + ^ > 0. 


(5.17a) 

(5.17b) 


It is readily seen that (5.17a I is satisfied for 9 > 1/4. For strictly positive £2 the inequality is 
satisfied whenever 9 > 1/4. Regarding inequality ( |5.17b ) we consider the left-hand side as a 
second-order polynomial in 9 with discriminant 

A = (2 + ^2? - 4(1 + ^2? = -^ 2 ( 3^2 + 4). 

If Z 2 > 0, then A < 0 and the polynomial is strictly positive for all real numbers 9. If Z 2 = 0, the 
polynomial reduces to 29^ — 29 + 1/2 which reaches its minimum (zero) in 9 = 1/2. 


Let 0 > 1/4 and 9 ^ 1/2. Then, applying Proposition |5.2| in this region yields 


lim |i?| = lim 
h^O h-fO 


+ pz + 9zoz + (^ - 9)z‘^ 


< 1 , 


and thus 


|i?^| = ^ VioO. 

Since Nq is independent from h and |1/(1— |z)| < 1 one may conclude that 

|t7^| = O(h^) Vw > 0. 

Next, consider the case 0 = 1/2. Recall that the MCS scheme then reduces to the original CS 
scheme. If \crj\ = it follows that 


lim 

^-0 4^sin"^ 


P = 9{1 + Z2), 
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with Z 2 > 0 such that proposition 5.2 can be applied and 


assume \cr]\ < h ® with 52 < 1/2. An expansion similar to (5.7) yields 


i?" I = O {h^) for all w > 0. Now, 


log 


-1 


le^sin^’^ 


{p^ +pz + ezoz+{^ -9)z'^) 


= log(-02 + 20 - 1) 


29 


1 


4Asin2 ’I 


— |0iai cot ^ — 9^z'^2^ + ^ (—pcot | iai cot 


— ^^ — \9iai cot ^ — 9 ^z2 ^ + 9pcoi 

— 2(1 — 0) (—pcot ^p + 1 ioi cot 

1 , I 0 I 


-02 + 20 - 


+ O 


sm 


2 i3i 


sm 


+ 0 ' 



and 


log 


le^sin^f 


= log 02 
+ 20 ( 
+ O 


4Asin2^ 2 


— ^0iai cot ^ — 02 ^ 2 ^^ I ^ 


02 


sm' 


2 'di 


+ 7] 


2 



Making use of ^ = 1 /2 it follows that 
^ I Jh' 


log(-i?) = 


Asin2^ 


2 ^ ' ^2 

-1 


h + O 


sin2 ^ I sin ^ 


+ p 


,2 



A sin2 ^ 


— Xp + \a 2 Xp \ h + O 


I 0 I 


sm2 ^ I sin ^ 


+ 0' 



and hence 


‘“dl-s)") = + “«) 


I 0 I 


sm - 


+ 02 /1 



In order to analyse the Rannacher time stepping we note 

1 , I0I 


1 — hz 

2 =1 + 0 


2Asin2’^ 


sin2 ^ I sin ^ 


+ 02 h 


which yields 


log (1 - = 2iVo log ^ 


I0I 


sin2^ 


Isinf 


+ 0 
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By exploring relationship (5.2) one may conclude that 


Un = (-1)^-^“ 


f^2No 


-1 


(2Asin"f 


exp 


A2 sin2 ^ 


- r] + ia 2 T] 


)(i+o((A 


JyL 

|,9i| 


(5.18) 

Whenever \cr]\ = h the right-hand side of (5.18) is O (h’") for all w > 0 such that we can use 
expression (5.18) for the whole region in the case of 0 = 1/2. 


5.6 Region 5: |k| < h |c? 7 | > h with qi < 1/2, q 2 > 1/2 


The analysis for this region is completely analogous to the analysis in Subsection 5.5 Hence, for 
0 > 1/4 and 0 1/2 it follows that 


\Un\ = 0{h^) Vw > 0 . 

Whenever the CS scheme is considered, i.e. 0 = 1/2, one gets the expression 


Un = 




(2A sin^ 


exp 


A2 sin2 ’I 


— K + \a2K 


(l + o((. 


l«l 

l>?2| 


(5.19) 


5.7 Connection with stability of the MCS scheme 

In the above analysis natural bounds on the MCS parameter 0 arise under which the asymptotic 
results are valid. These bounds can be interpreted as stability bounds. In particular, the condi¬ 
tions 0 > 0 > are needed to ensure that the Fourier transform Un is negligible in the 

second region. This restriction is only slightly stronger than the lower bound on 0 derived in [^, 
guaranteeing unconditional stability of the MCS scheme in the von Neumann sense pertinent to 
two-dimensional diffusion equations with mixed derivative term. This is, indeed, not very sur¬ 
prising. In [5] it is stated that the stability analysis of the MCS scheme in this case reduces to 
bounding by one of the modulus of the scalar expression 

z {0zo + - 0)T)Z 


where z = zq + zi + £ 2 , p = {1 — 0zi)(l — 0 Z 2 ) and zq, £[, £2 denote real numbers satisfying the 
condition ( 5.12[ ). This explains why Proposition 5.1 is just a slight modification of one of the 
statements in 0 Theorem 1]. 


6 Asymptotic analysis in physical space 

In this section we will use the asymptotic results in Fourier space from Section to perform an 
error analysis in physical space. First note that the Fourier transform u is only sizeable in region 
1 of the Fourier domain. In the other regions it holds that k > or crj > and hence 


u{K,r],l) = O (h^) Ww > 0. 
Based on equalities (|3.2[), (5.10) and (|5.16|) we define 


Etow _ ^2 exp (—^2 _ 2 pK'q — r]'^ + iaiK + m2?7)(s^^^ (k, rj) + Nq^{k, r])) 


r[2]/ 


and 


^high 


( 1 t(l?l,-02) \ 
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Recall that i?i = nhi, '§2 = * 7/12 and /i 2 = chi = ch. As a consequence, E^°'^ is only sizeable in 
region 1 and is only sizeable in region 3. In the other regions of the Fourier domain Ujsi is 

negligible whenever 9 > max{|, and 6 ^ 1/2. Hence, for these values of 9, the results can 

be combined to 

77 / 12 ) — m(k,? 7 ,1) « |k|, |c77| < 7r//i. (6.1) 

When 9 = 1/2, i.e. when the MCS scheme reduces to the CS scheme, Uiq is also sizeable in region 
4 and region 5 of the Fourier domain. This case will be treated separately. 


6.1 MCS scheme with 6 ^ 1/2 

Consider the case where the MCS scheme is different from the CS sche me, i.e. 9 ^ 1/2, and 
suppose that the restriction 9 > max{i, is satisfied. Approximation (6.1 1 is then valid and 

based on (5.41 we have for the total error: 


where 


fc - w(a7„ 7/fc, 1) « E\°^ + Eft. 

l2 /*oo noo 

= ^ / / u(K,77,l)(s[2l(7c,77) + 7V^^'(K,77))exp(iraj)exp(i77yfe)dKd?7 (6.2) 

J —00 J —00 


and 


E 


high ^ fe 2 JVo^ 4 JVo |■K/h2 n^/hi exp(iKa;j) exp(i 777 /fc) 

4'^^ ■n-//i2 tt/Tii [2A7.(75i,' d2)] ° 

l^2No-2^ANo-i /-TT exp(ij'di) exp(iA:-d2) 


j,k 


exp - 


t(7?i,'d2 


4A20^ sin^ t sin^ f- 


47r^ 


I -TT J -TT [2At(7?l, 7 ^ 2 )] 


2No 


exp 


1 




4A202 sin 2 4 sin 2 


dKdrj 


d'didd2- 


First, consider the low-wavenumber error. The inverse mixed discrete/continuous Fourier trans¬ 
form of is given by 

^7r//i2 pT^lhi 

/ / u{K,ri,l){s^‘^^{K,ri) + Nf{K,r]))exp{iKXj)exp{ir]yk)dKdr], 

'—TT{h2 J—Trlhi 


47r2 


which can be approximated by (6.2) as hi,h 2 tend to zero. Let 

/ 


(l>pix,y) = 


•y/47r2(l-p2) 


exp I 


— 2pxy-\-y^ 

2(1-P") 


the density function of a two-dimensional standard-normally distributed random variable with 
correlation p. Its Fourier transform is 


^p{k,v) = exp - pKT]-\ 




Hence, for all positive integers ni,n 2 the Fourier transform of 


gni+>i2 , / x+ai y+a2 ^ 

dx'-idy'^^Vp y V2 ’ \/2 / 

2 (i-\/ 2 K)"^ (i-\/277)"^ exp(—— 2pKp — + iaiK + ia 2 p), 

such that the inverse Fourier transform of 

exp(—— 2pKp — + ioiAt -I- \a 2 p)K"^p"^ 

is given by 


IS 


1 


gni+n2 


2 (iy2)"i+"2 dx^^dy” 


X + ai y + a 2 
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Recalling the formulas for and from Subsection 5.2 this leads to the following expression 
for the low-wavenumber error: 

(6.3) 


rplow _ p.2/^low 


with 


C‘. 


low 

^j.Vk 


1 P f 

48 9a:'^ 12 \dx^dy 




-b 


1 52 

2 dx'^ 

A2 /I 

2 dx"^ 

1 


12 

A2 


oi d 
\f2 dx 

^ dxdy 

92 


2 dx"^ ^ dxdy 

92 


2 i9a;9y 

4 12 


254 


1 


02 ^ \ 

V2dy) 


dxdy^) 

1 92 

2 

1 ^2 

2 9j/2 ^2 

1 ^2 

2 9j/2 ^2 

1 92 


c 2 94 
48^ 


«! 


^3 


I2V2 5a:3 12^2 dy^ 


020^ 






1 52 

2 2 c}?/2 


iVnA2 


52 

^ dxdv 


1 92 

2 


Cl d 
, /o 


1 52 

2 9a;^ 

02 5 y 

V2dy) 
a2 d\ 
V2dy) ^ 


oi 9 
y/2 dx 

02 5 y 


92 

dxdy 


1 92 

2 dy'^ 


ai d 02 9 
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Next, consider the high-wavenumber error and note that 

exp(iji?i) exp(iA:d 2 ) = cos(ji?i-b fci? 2 ) + i sin(ji?i-b fct? 2 )- 


Symmetry yields 
where 
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Combining both expressions (6.3) and (6.4) gives an approximation for the total error: 


UN,j,k - uix„yk, 1) « 


(6.4) 


(6.5) 


The values are only dependent on the position {xj,yk) = {jhi,kh 2 ), the parameter values 

of the problem and the ratios c and A. The constants only depend on the index {j,k), 

the correlation parameter p and the ratios c, A. For the numerical experiments, cf. infra, the 
values are calculated by determining all the partial derivatives. The integrals in are 

approximated by numerical integration. It is readily seen that 

maxiyf | = |<ri. 


so has a maximum magnitude where {xj,yk) = (0,0). This is exactly at the position of 

the discontinuity of the initial function. At the end of Subsection |5.4| it was conjectured that 
for larger values of the MCS parameter 6 one can expect a larger high-wavenumber error. This 
conjecture is confirmed by the above analysis given that (.('^ 1 ,^ 2 ) is always positive. In order 
to avoid spurious erratic behaviour in the numerical solution, it is therefore recommended to use 
smaller values of the parameter 0. However, one has to take into account the lower bound on 0 
described in Subsection Ell 
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We showed in (6.5) that the total error is jVo = 2 is a lower bound 

on Nq for the Rannacher time stepping in order to ensure convergence of the numerical solution 
to the exact solution. This is confirmed by the plots in Figure which display total errors (in 
the maximum norm) in actual numerical experiments for model problem (3.1) as a function of 


l/h, with parameter values p = —0.7, oi =2,02 = 3, MCS parameter 9 = 1/3 and with c = 1, 
0.2 < A < 0.8. Since it is not possible to handle infinite domains in numerical experiments, the 
computational domain is restricted to spatial gridpoints {xj,yk) G [—10,10] x [—10,10]. At the 
boundaries, homogeneous Dirichlet boundary conditions are applied. In the left plots the case 
Aq = 0 is considered, whereas the right plots show the corresponding results for Nq = 2. In the 
upper plots the maximum error between our numerical solution and the exact solution is shown as 
a function of l/h for different values of A. In the lower plots we show the same maximum error for 
one value of A, together with our theoretical estimates for the corresponding low-wavenumber error 
and high-wavenumber error. In these lower plots it is clearly seen that our theoretical estimates 
for the total error are sharp. 

For the case where no Rannacher time stepping is applied, the left plots in Figure [^ reveal 
second-order convergence behaviour until h reaches a critical value where the high-wavenumber 
error starts exceeding the low-wavenumber error. It can be observed that this value of h, and thus 
the high-wavenumber error, is highly dependent on the ratio A = At/h. For smaller values of A, 


E- 


high 


f, is only sizeable whenever h is very small, whereas for larger values of A, ^ already dom¬ 
inates the total error for larger values of h. Moreover, the error constant for the low-wavenumber 
error is also dependent on A. However, this is much less pronounced than for the high-wavenumber 
error. 

The right plots in Figure show the corresponding results in the case where the first two MCS 
timesteps are replaced by four backward Euler half-timesteps, thus Nq = 2. One observes that 
the numerical approximations now exhibit second-order convergence for all values of A. In the 
bottom right plot the high-wavenumber error is not visible since it is strongly dominated by the 
low-wavenumber error. The same observation is made for other values of A. Hence, whenever 
Rannacher time stepping is applied with Nq = 2 , the total error can be approximated by 
which is of second-order in h. We find that the error constant for the low-wavenumber error is 
mildly dependent on the ratio A = At/h. This can be explained through the fact that for a fixed 
value of h but smaller value of A the same semidiscrete system is solved with a smaller timestep 
At. Finally, we notice that the latter error constant is slightly larger than for the case where 
Nq = 0. Thus, by applying Rannacher time stepping with Nq = 2, second-order convergence can 
be recovered at the small cost of a marginally larger error constant for the low-wavenumber error. 

As stated above, the high-wavenumber error is very sensitive to the MCS parameter 9. To 
illustrate this, Figure]^ shows the same plots as in Figure]^ but with the MCS parameter replaced 
by 0 = 1. It can be seen that all the conclusions from Figure [^remain valid. In order to get decent 
plots, however, it is necessary to consider smaller values for A. This confirms that, for fixed A, 
is strongly increasing as a function of 9. Moreover, by comparing the upper plots of Figure 
[^and Figure]^ for A = 0.2 it can be seen that the error constant of the low-wavenumber error is 
substantially larger for MCS parameter 9 = 1 than for 9 = 1/3. We conjecture that for fixed A 
and fixed h, the low-wavenumber error is also increasing as a function of 9. Therefore, regardless 
of the number of Rannacher timesteps Nq , it seems more favourable to consider smaller values of 
9. In particular, the lowest value of 9 which satisfies the restrictions from Subsection |5.7| for all 
values 1 p 1 < 1 is given by 0 = 1/3. 


6.2 MCS scheme with 9 = 1/2 

For 0 = 1/2, the MCS scheme reduces to the CS scheme and Un is not negligible in region 4 and 
region 5. Based on equalities (5.18) and (5.19) we define 
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Figure 5: Convergence of the numerical solution for Nq = 0 (left) and Nq = 2 (right). The 
parameter values are: p = —0.7, Ui = 2, 02 = 3, 0 = 1/3. 
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Since di = Khi, ^2 = 77/12 and /12 = chi = ch, Qjjiy considered in region 4 and 

j^CS,b jg Qjjjy jjQ^ negligible in region 5 . Hence, the Fourier error ( 5 . 3 ) can be approximated by 

UN{Khi,'qh2) - u{k, 77,1) « |k|, |c?7| < ir/h. 

The inverse mixed discrete/continuous Fourier transform of + is given by 
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As hi , /i2 tend to zero this can be approximated by 
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Making use of a symmetry argument and a one-dimensional inverse Fourier transformation, Ef'^ 
can be rewritten as 
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where (j) denotes the density function of a standard normally distributed random variable. It is 
readily seen that respectively reaches its highest magnitude near the points {j,k) 

where {xj,yk) « (0,-02), respectively {xj,yk) ~ (—ai, 0 ). For the numerical experiments, the 
integrals in Cj;® and are approximated by numerical integration. Combining the expressions 
( 6 . 3 ), ( 6 . 4 ) and ( |6.6[ ) leads to the following approximation of the total error: 

(6.7) 




From approximation ( 6 . 7 ) it can be concluded that the total error is also O ( 7 i™'"{ 2 , 27 Vo- 2 }^ 
when CS time stepping is considered. This matches the observations from the plots in Figure 
which show convergence results for the same problem as in Subsection |6.1| but with MCS parameter 
6 = 1 / 2 . The lower plots indicate again that our theoretical estimates for the total error are 
sharp. Without Rannacher time stepping, i.e. Nq = 0 , the results in Figure show second-order 
convergence in h until E^^ starts exceeding the low-wavenumber error. Then the total error 
increases in a first order way until the high-wavenumber error starts dominating. From there the 
total error is O {h~^). In case the MCS scheme is replaced in the first two timesteps by four 
half-timesteps of the implicit Euler scheme, i.e. Nq = 2 , Figure reveals unconditional second- 
order convergence in h. Note that both E^^ and are not visible in the lower-right plot 

because they are strongly dominated by the low-wavenumber error. The same observation as in 
Subsection | 6 . 1 | can be made concerning the dependency of the low- and high-wavenumber error 
on the parameter A. 
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7 Conclusion 


If the initial data is nonsmooth, application of the MCS scheme for multidimensional time- 
dependent convection-diffusion equations with mixed-derivative terms can cause spurious erratic 
behaviour in the numerical solution. A motivating example, with the two-dimensional Black- 
Scholes equation for a two-asset cash-or-nothing option, shows that this undesirable feature can 
be resolved by replacing the very hrst Nq MCS timesteps by 2No half-timesteps of the implicit 
Euler scheme, with Nq = 2. We proved, by Fourier analysis, that for a model two-dimensional 
convection-diffusion equation with mixed-derivative term and with Dirac delta initial data, the 
total error can be approximated by the sum of a low-wavenumber error of 0{h?) and a high- 
wavenumber error of In case the MCS scheme reduces to the CS scheme, i.e. when 

6 = 1/2, this has to be augmented with an extra error term of Hence, Aq = 2 is the 

minimum on Nq in order to guarantee (second-order) convergence of the numerical solution to the 
exact solution, in the maximum norm. In general this choice for Nq is optimal since larger values 
will increase the low-wavenumber error. Our convergence analysis and numerical experiments fur¬ 
ther indicate that it is favourable to consider small values of the MCS parameter 9. However, it 
is necessary to take into account the lower bounds on 6 in order for our asymptotic analysis to be 
valid. The smallest value which satisfies all the restrictions, independent of the parameters of the 
model, is given by 0 = 1/3. This is, indeed, also the most common value of 9 for the MCS scheme 
considered in the literature. 
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